Hierarchical Liouville-space approach for accurate and universal characterization of 

quantum impurity systems 
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A hierarchical equations of motion based numerical approach is developed for accurate and efficient 
evaluation of dynamical observables of strongly correlated quantum impurity systems. This approach 
is capable of describing quantitatively Kondo resonance and Fermi liquid characteristics, achieving 
the accuracy of latest high-level numerical renormalization group approach, as demonstrated on 
single-impurity Anderson model systems. Its application to a two-impurity model reproduces the 
experimentally observed crossover behavior, which affirms the continuous transition from the Kondo 
singlet state of individual impurity to the singlet spin-state formed between two impurities. 
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Quantum impurity systems cover a broad range of im- 
portant physical systems where strong electron-electron 
(e-e) interactions among a few localized impurities affect 
crucially the system properties. Besides the e-e interac- 
tions, the impurities are also coupled to the itinerant elec- 
trons in surrounding bulk materials, which serve as the 
electron reservoir as well as thermal bath for the impu- 
rities. Moreover, some extensive strongly correlated sys- 
tems can be treated as quantum impurity systems. For 
instance, the celebrated Hubbard model can be mapped 
onto an Anderson impurity system via a self-consistent 
dynamical mean- field theory The strong e-e inter- 
actions give rise to a variety of intriguing phenomena 
of prominent many-body nature, such as Kondo effects, 
Mott metal-insulator transition, and high-temperature 
superconductivity. Examples of localized impurities are 
the d- or /-electrons of transition metal atoms, and elec- 
trons trapped in quantum dots. 

Accurate characterization of quantum impurity sys- 
tems is the key to understand the mechanisms and ef- 
fects of strong electron correlations. This has remained 
a very challenging task, especially for the quantitative 
evaluation of dynamical quantities which are directly re- 
lated to experimental measurements, such as the pro- 
jected density of states and spectral function of the local- 
ized impurities. A vast amount of theoretical efforts have 
been devoted to achieving this goal. A variety of numeri- 
cal approaches have been developed, including the quan- 
tum Monte Carlo approach 0] , density matrix renormal- 
ization group method [3L numerical renormalization 
group (NRG) method p, |6(, many-body perturbation 
theory Q, and effective/quasi single-particle approaches 
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[1, Q . Despite their success in elucidating some funda- 
mental features of strongly correlated systems, the prac- 
ticality of existin g app roaches has been limited within a 
few basic models [1 0l — ll 2j j . The reason is mainly twofold: 
(?) the applicability of involving techniques relies criti- 
cally on the system configuration, and (ii) the complexity 
of numerical algorithms employed increases dramatically 
with the number of impurities. Consequently, generaliza- 
tion of existing approaches to more complex models 
is often very difficult. Therefore, an accurate and uni- 
versal approach capable of addressing strong correlation 
effects in general quantum impurity systems is highly de- 
sirable. 

In this Letter we propose a general approach based on 
a hierarchical equations of motion (HEOM) formalism 
[lH to characterize quantum impurity systems from the 
perspective of open dissipative dynamics. Using this ap- 
proach, a variety of dynamical observables of significant 
experimental relevance can be evaluated accurately and 
efficiently. We first exemplify its practicality through the 
evaluation of spectral functions of single-impurity Ander- 
son model (SIAM) systems. Its accuracy is shown to be 
comparable to the latest high-level NRG approach. The 
associated Kondo resonance and Fermi liquid character- 
istics are scrutinized by comparing directly to existing 
numerical and analytical results. Calculations are then 
extended to a two-impurity Anderson model (TIAM) sys- 
tem, where unconventional Kondo resonance features is 
disclosed which well reproduces a recent experiment [14J . 

In the framework of HEOM, the localized impurities 
constitute the open system of primary interest, while the 
surrounding reservoirs of itinerant electrons are treated 
as environment. The total Hamiltonian consists of the 
interacting impurities (H sys ), the noninteracting elec- 
tron reservoirs (H rQS ), and their coupling which assumes 
the form of i? sys - res = Z^fcCWfc d^d ak + H.c). Here 
d' u and da denote the creation and annihilation opera- 
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tors for impurity state \fi) (including spin, space, etc.), 
while d/, and d a k are those for the a-reservoir state \k) 
of energy e Q fc. The influence of electron reservoirs on 
the impurities is taken into account through the overall 
reservoir spectral functions, A^(u>) = J2 a A QA1 „(u;) = 
"■Eaic WCit % ~ £ q/c), in the absence of applied 
chemical potentials. The HEOM that govern the dynam- 
ics of open system assume the form of [1 31 ] : 
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where the basic variables are the reduced system density 
operator p^(t) = tr rcs /Ototai(i) and auxiliary density op- 
erators (ADO), P^..j n (t); n = 1, • • • , L, with L denoting 
the terminal or truncated tier level. The Liouvillian of 
impurities, £■ = h^ 1 [H sys , •], may contain both e-e in- 
teraction and time-dependent external fields. The super- 
operators Aj and Cj are expressed by Eqs. (SI) — (S2) of 
Ref. [Hj]. Here, the individual index j = (er/im) corre- 
sponds to electron transfer to/from (a = +/—) the im- 
purity state fi associated with a memory time 7~ . The 
total number of distinct j-indexes involved is determined 
by the preset level of accuracy for decomposing reservoir 
correlation functions by exponential functions. Such a 
number draws the maximum tier level L max at which the 
HEOM of Eq. (fTJ ultimately terminate pf . 

The HEOM formalism is nonperturbative. It has been 
demonstrated that the hierarchy is self-contained at L = 
2 level for noninteracting H sys [l3[ ■ In practice, Eq. (fT]) 
needs to be truncated at a certain L, followed by a conver- 
gence test. Extensive numerical tests have affirmed that 
a relatively low truncation tier (L 4) is usually suffi- 
cient to yield quantitatively converged results for weak 
and medium impurity-reservoir couplings. The HEOM 
approach has been applied to study the static and dy- 
namical electron transport properties of quantum dot 
systems, with which some interesting phenomena have 
been revealed, such as the dynamical Coulomb blockade 
and dynamical Kondo current (l7j . 
For evaluation of dynamical variables of the impuri- 
ties system, we focus on correlation function between two 
arbitrary dynamical operators, Cab(£) = (A(t)B(0)} = 
tatotai[^(i)-B(0)p^ tal (T)]. Here, the Heisenberg opera- 
tors and thermal equilibrium density operator /? total (-^) 
are all defined in the total space. A linear response the- 
ory in the HEOM space is established (see Supplemen- 
tal Materials [l5[ for details), based on which Cab {t) 
is retrieved exactly within the HEOM framework. Let 
CUb(w) = i J dte lu>t CAB(t), which satisfies the de- 
tailed balance relation of Cba(—u) = e~ hu) / kBT Cab(w). 
The system spectral function is obtained as Jab{u) = 
l 
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Jdte^({A(t),B(0)}) 
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= l(l + e-^/ k - T )C A B(Lo). 
it recovers the spectral func- 
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FIG. 1: (Color online). The spin-up or down spectral function 
of an asymmetric SIAM calculated by the HEOM approach 
at different truncation tiers. The inset magnifies the Kondo 
resonance peak at ui = 0. The parameters adopted are ed = 
-5, U = 15, W = 10, and T = 0.075 (in unit of A). 



tion of the impurity state /i, i.e., A^(uj) = J a -t (w) = 

— — ImG^^a;), which can be experimentally measured 
via angle-resolved photoemission spectroscopy [l8| and 
scanning tunneling microscope (ljjj . Here, G MM (w) is the 
retarded Green's function. 

For numerical demonstrations, consider first an asym- 
metric (U ^ —2 e d ) SIAM system of H sys = £d("-t + ?V) + 
Un^-fii, where n M = ct^a^. The reservoir spectral func- 
tion A MJ/ (w) = S^AW 2 j '(uj 2 +W 2 ) is adopted, where A is 
the effective impurity-reservoir coupling strength and W 
the reservoir band width. Such a model has been widely 
studied by various methods [20j, including perturbation 
theory and NRG approach. Figure [T] shows the calculated 
impurity spectral function A(u>) via the HEOM approach 
up to the converged tier level. The parameters adopted 
are (in unit of A): ed = —5, U = 15, W = 10, and 
T = 0.075. The calculation results display well-known 
features of SIAM: (i) The two resonance peaks resolved 
at around u) — ed and U + ed represent the single-electron 
and two-electron occupancy states, respectively; in) The 
peak at the Fermi energy (uj = Ep = 0) highlights the 
presence of Kondo resonance under a low temperature; 
(Hi) The sum rule J A(uj) du — 1 is satisfied to numeri- 
cal precision. The comparison in Fig.[T]dcmonstrates dis- 
tinctly that the HEOM results converge rapidly with L 
for full energy range. This confirms the HEOM approach 
is capable of yielding quantitatively accurate dynamical 
properties at relatively small cost. 

Figure [5] depicts the calculated A(lu) of a symmetric 
(U = -2e d ) SIAM with T = 0.2A and W = 50A, from 
weak (U = 0.57rA) to strong (U — 6irA) e-e interaction. 
For comparison, we also show results obtained by using 
the full density matrix NRG method [2l|, where a self- 
energy scheme of Ref. [22( is employed, and the results 
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FIG. 2: (Color online). Spectral function of a symmetric 
SIAM calculated by HEOM and NRG at different e-e inter- 
actions. The inset shows the imaginary part of interaction 
self-energy calculated from HEOM at energy close to to — 0. 
Note the break of y axis in the inset. 



are averaged over 8 different logarithmic discretizations 
[Hj]. Note that our NRG curves in Fig.H differ slightly 
from those in Ref. [24| , due to the difference in lineshape 
of A^i/(a;) (Lorentzian versus constant). Apparently, the 
HEOM results agree quantitatively with the NRG coun- 
terparts for all e-e interaction strengths studied. In the 
weak (U = 0.5 and 1.0 7rA) and strong (U — 67rA) inter- 
action regimes, HEOM and NRG curves almost overlap 
with each other. Minor deviation is observed at an in- 
termediate value (U — Sir A), where the Hubbard peaks 
(around to = and + U) from HEOM are at slightly 
higher energy than those from NRG. As is well known, 
NRG method tends to push the Hubbard peaks towards 
the lower energy, due to the Hilbert space truncation and 
log-Gaussian broadening. Therefore, the observed devi- 
ation takes place in the expected direction. The small 
differences near u> = for intermediate and strong in- 
teractions may result from the different numerical tech- 
niques adopted in NRG and HEOM approaches. 

Highlighted in the inset of Fig. [5] are the imaginary part 
of interaction self-energy (scattered circles) , exhibiting a 
parabolic lineshape near lu = Ef = (solid lines). This 
is a clear indication of Fermi liquid character [2(| ■ It has 
been proved by Luttinger [25[ that for a symmetric SIAM 
at T = 0, the height of Kondo-resonance peak is exactly 
l/VA, independent of U. In general, at finite T and U, it 
is expected that A(ui = 0) < l/VA for both asymmetric 
and symmetric SIAM systems |26|, as affirmed in Fig.Q] 
and Fig.[2j respectively. 

We extend the generic HEOM calculations to a 
parallel-coupled TIAM system, where the two-impurity 
Kondo effect emerges. We have now H sys = Hi + H 2 + 
V12, with Hi (H2) being the SIAM Hamiltonian for the 
impurity 1 (2) and V\<z — t(a{*a2t+ a ii a 24.+H.c). Such a 
TIAM model has been realized experimentally via a dou- 
ble quantum dot system as sketched in FigJ^a), where 
the inter-dot coupling t is tuned by plunger gates [l4| . 
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FIG. 3: (Color online), (a) Schematic diagram of a parallel- 
coupled TIAM system, (b) Spectral function of a symmetric 
TIAM consisting of two identical impurities, A(u) = Ai(uj) — 
A2(lo), at various inter-impurity coupling strengths t ranging 
from to 2.5A. (c) HEOM results and (d) experimental mea- 
surement of A{lu) in Ref. [3 (Reprinted with permission), 
where the red arrows indicate the increasing direction of t. 



An intriguing question is how the inter-impurity cou- 
pling would affect the Kondo signature. It is known 
that the i-coupling should induce a transition from a 
Kondo singlet of individual impurity to a singlet spin- 
state between two impurities. However, it remains un- 
clear whether such a transition involves a critical point 
(quantum phase transition) or not. Jones et al. [27j stud- 
ied a serial-coupled two-impurity system by carrying out 
pioneering NRG calculations on a Kondo model with an 
impurity-impurity exchange coupling term. Their results 
marked the transition by a non-Fermi liquid fixed (crit- 
ical) point. Recent NRG calculations also found the 
critical point in a parallel-coupled two- impurity Kondo 
model 28]. However, such a quantum critical behavior 
is yet to be verified by experiments. In contrast, Chen 
et al. have observed a continuous crossover behavior in a 
parallel-coupled double quantum dot setup represented 
by Fig.[3][a) [IH . The apparent inconsistency between 
theory and experiment calls for quantitatively accurate 
numerical studies on more realistic models. For this pur- 
pose, we explore a TIAM by carrying out HEOM calcu- 
lations to numerical convergence. 

Figure [3jb) summarizes the calculated spectral func- 
tions of a TIAM consisting of two identical impurities, 
A(uj) = Ai(uj) — A2(uj), at various values of t. Other 
parameters are (in unit of A): W — 2.5, U% — U2 = 10, 
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€\ = e 2 = —5, and T = 0.5. For a symmetric TIAM, the 
differential conductance (dl/dV) is approximately pro- 
portional to A(ui) in linear response regime at low tem- 
perature [lij]. Therefore, the HEOM calculated A(u>) can 
be compared directly to experimentally measured dl /dV 
of Ref. [14} : see Fig.[3tc) and (d), where red arrows indi- 
cate the increasing t. With Fig.[3Jb)-(d), we demonstrate 
that the calculations well reproduce the essential exper- 
imental features: (i) The system undergoes a transition 
from a Kondo singlet involving individual impurity (char- 
acterized by the single-peaked lineshape) at t < A, to the 
singlet spin-state between two impurities (characterized 
by the double-peaked lineshape) at t > 1.5 A; (ii) The 
transition exhibits a continuous crossover. As t increases 
from A, the single Kondo peak first broadens and ap- 
proaches to its maximal height before it drops and splits 
into two; (Hi) Beyond the crossover regime (t > 1.5 A), 
the separation between the two peaks increases rapidly 
with t, with each peak resembling the original unified one 
in both shape and width. 

For a TIAM system, the nonzero t gives rise to an ef- 
fective anti-ferromagnetic coupling, J = At 2 /U , between 
the local spin moments at the two impurities. At a weak 
J, the two spin moments are nearly independent of each 
other, and the itinerant electrons screen the local spin 
at each impurity separately. In contrast, at a sufficiently 
strong J, singlet spin-states are formed, which span over 
both impurities. The absence of quantum phase transi- 
tion in Fig. [3] might be due to the instability of quantum 
critical point under the chosen parameters, or caused by 
the charge transfer term Vn in our model. The latter 
may change the possible quantum phase transition to a 
continuous crossover [29|, [3Cj . More comprehensive inves- 
tigation is needed. 

To summarize, the practicality of our developed hierar- 



chical Liouville-space approach is demonstrated through 
studies on Anderson impurity models, where the key 
Kondo resonance and Fermi liquid features due to strong 
e-e interaction are accurately characterized. The HEOM 
approach can be straightforwardly extended to more gen- 
eral quantum impurity models without additional deriva- 
tion and programming efforts. Besides the equilibrium 
dynamical observables, the HEOM approach is also ca- 
pable of giving transient electronic response to external 
fields such as laser pulses or applied voltages 0, Il7} . 
These make the HEOM approach a powerful tool for 
simulation of strongly correlated quantum impurity sys- 
tems. We should point out that the HEOM calculations 
can be expensive at extremely low temperature or strong 
impurity-reservoir coupling, where a high truncation tier 
level is necessary to achieve quantitative convergence. It 
is however possible to reduce computational cost signif- 
icantly by developing more suitable reservoir spectrum 
decomposition schemes. Once converged, the HEOM re- 
sults can serve as benchmarks to calibrate various ap- 
proximate but efficient numerical approaches, particu- 
larly the effective single-electron approaches, which are 
useful for studying more complex systems. Moreover, it is 
anticipated that HEOM would become a routine numer- 
ical solver for strongly correlated lattice systems within 
the framework of dynamical mean field theory [l[. 
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